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1. Introduction 

The existence of singular points at which the partition function (JT) vanishes is expected in 
a complex parameter plane. These are called Lee- Yang zeros or Fisher zeros |jl|]. The scaling 
analysis of such singularities in the complex parameter plane is an interesting approach to under- 
stand the nature of phase transitions. A few years ago, the thermodynamic singularities of QCD 
in the complex chemical potential (ju 9 ) plane were discussed by M. Stephanov using universality 
arguments in the vicinity of the critical point and a random matrix model It is interesting to 
perform numerical simulations of lattice QCD and to compare the results with the predictions from 
the universality arguments. 

Moreover, the study of the singularities in the complex plane is important to estimate the radius 
of convergence of a Taylor expansion in terms of \i q . The Taylor expansion method is widely 
used for the study of the equation of state at finite \i q . The expansion coefficients of pressure 
are defined by p(T ,H q ) /T = \n£>f /(VT 3 ) = Y,n c n{T){pL q /T) n , where V is spatial volume and 
T is temperature. Because p is an analytic function of \i q for finite volume, c n does not change 
even when \i q is a complex number. Hence, the convergence radius is determined by the nearest 
singularity from \i q = in the complex \i q plane. 

In this report, we study the thermodynamic singularities of QCD by a simulation with rela- 
tively heavy quark mass. Although simulations near the chiral limit are required for the study of 
universality class in 0, it is important to study the distribution of the singularities in the complex 
jj, q plane as a first step toward the universality argument and for an estimation of the convergence 
radius even if the quark mass is not very small. In Sec. ||, we discuss a method to investigate the 
singularities at complex jj, q and point out problems in this method. The properties of Lee- Yang 
zeros are also discussed, introducing the probability distribution function of complex phase of the 
Boltzmann weight. In Sec. ^ we estimate the distribution of the singular points using the data 
obtained by a simulation of two-flavor QCD with p4-improved staggered quarks of m % lm p « 0.7 
on a 16 3 x 4 lattice in [[J, and discuss the convergence radius. Conclusions are given in Sec. [J 

2. Singularities and complex phase of quark determinant with complex \i q 

The grand partition function at a complex chemical potential \i q = jitR e + ipi\ m is given by 

^(j3(r),^ Re + ^ Im ) = J @U(detM(p: Re + ip: lm )) N! e- s s. (2.1) 
To find the singularities in the complex plane, we define a normalized partition function, 



•^norm (P , £*Re + « Hba ) = 



^(j8,jURe+z'Atln 



( detM^Re + Z^in 

V detAf(O) 



(J3,M„=0) 



(2.2) 



where (••■)(/? /i»=0) means the expectation value at \i q = 0, and j3 = 6/g . Because 3f(j5,0) is 
always nonzero for finite volume, the position of Jf(j8,jU 9 ) = can be found by calculating £¥ norm . 
Since £t? vanishes due to the complex phase of detM, i.e. 6 = A^fIm(lndetM(/i^)), we introduce 
the complex phase distribution function W (8). The partition function is then written by 



•^lorm (fi i MRe + ' Mln 



e w W(d)dd 



(2.3) 



2 



Singularities of QCD in the complex chemical potential plane 



Shinji Ejiri 



Identifying 6 + 2nn with 6, 3f novm vanishes when the contributions from 6 and 6 + % cancel each 
other. Using the distribution function, we point out that the phase 8 contains two components: one 
is related to the total quark number and produces characteristic properties of Lee- Yang zeros, and 
the other is irrelevant to Lee- Yang zeros and causes the sign problem at large \l q . 

We discuss the complex phase in the vicinity of the real ji q axis. The phase is given by 



0(jURe + /jUlm) =Nf 



(2.4) 



Because 3f is real and positive for real \i q , the \i\ m -independent term does not contribute to 3? = 
for finite V. However, because of statistical fluctuations in the /ii m -independent part of 6, 3f 
may be smaller than the statistical error even at pL\ m = 0. We plot the histogram of the complex 
phase in Fig. [I] (left). The phase is calculated by a Taylor expansion of IndetM up to 0{ji q ) using 
data in a simulation with p4-improved staggered quarks [St]. The result of pL\ m /T = is given 
by the /x Im -independent part only. The magnitude of 3f decreases exponentially as ;U Re increases, 
i.e. 3f ~ exp[-(0 2 )/2] with (d 2 ) ~ 0(^| e ) for jj, lm = [§. Once 3f becomes smaller than the 
statistical error at large jl q , 3f n0Tm vanishes at random. Such 3f n0Tm = are irrelevant to Lee- Yang 
zeros, but one cannot distinguish such fake zeros from real ones. This is a kind of the sign problem. 

Another interesting point is that the operator d(]ndetM(n q ))/d(n q /T) in the second term of 
Eq. ( |2.4[ ) corresponds to the quark number (N) on each configuration. Therefore, the distribution 
of the complex phase is related to the distribution of the quark number at the leading order of Hi m . 
Let us consider a canonical ensemble with fixed N. The canonical partition function 3?c(T,N) for 
each N is related to the grand partition function 3f(T,pL q ) through a fugacity expansion, 

3f(T,^ q ) = £ 3f c {T,N)e N ^ /T = £ W(T,N)e i§ = f VW{T,Vp)e ipv ^ T dp, (2.5) 

where p is the quark number density. In the case of \L q = /^R e + i\l\ m , the complex phase in £¥(T, pL q ) 
is 6 = NpL\ m /T = pV/J.i m /T, since 3?c(T,N) is real and positive. W(T,N) is the distribution func- 
tion of the total quark number, W(T,N) = ^ c (T,N)e N ^/ T . 

The jU Im -independent part in this is eliminated if the canonical partition function is obtained 
by a partial path integral with fixed N. Once the problem of the unnecessary part of the phase is 
solved, one can discuss the order of phase transitions in the following way. 



We find from Eq. (|1J) that i2° as a function of (V}Xi m /T) is obtained through a Fourier trans- 
formation of W with respect to p . In the case of a normal point of jitR e or a crossover pseudo-critical 
point, where the distribution is expected to be of Gaussian for sufficiently large V, 3? does not van- 
ish except in the limit of VpL\ m /T — > ±oo because the function which is obtained through a Fourier 
transformation of a Gaussian function again is a Gaussian function. Therefore, the complex fi q at 
which 3? = does not go to the real axis in the infinite volume limit. 

On the other hand, at a first order phase transition point, two phases having a different quark 
number coexist. In this case, we expect that W(T,Vp) has two peaks having the same peak height at 
the transition point. Performing the Fourier transformation of such a double peaked function leads 
to a function which has zeros periodically. For example, a distribution function W(T, Vp) having 
two Gaussian peaks at p\ and P2 leads to 3f which has zeros at }l\ m /T = (2n + Y)n/\V{p2 — Pi)], 
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Figure 1: Left: The histograms of complex phase 9 for complex flq/T = HR e /T + ijii m /T at j3 = 3.65. The 
dashed curves are Gaussian functions. Right: The plaquette histogram w(P,P) and the effective potential 
Veff (P, J3 , 0) at \L q = for each j8 . 



with n = 0, 1,2,3, ■ ••. The Lee- Yang zeros approach the real axis as 1/V. Therefore, the study of 
1 /V scaling of the Lee- Yang zero is equivalent to finding the double-peak structure at a first order 
phase transition. The same discussion for SU(3) pure gauge theory is given in [g]. 

3. Grand canonical partition function by a density of state method 

As noted above, the numerical study in terms of the singularities of ^ = has a potential 
danger for large ji q . However, the position of = can be estimated from the distribution function 
of the complex phase, since 5° vanishes when the distribution function has two peaks or more and 
the contributions from these peaks cancel each other. In this study, we investigate the distribution 
function of the complex phase instead of iF. Once the quark number is fixed, the distribution 
function is related to the canonical partition function, as discussed in [||] for real \i q . However, 
fixing the quark number is not essential for eliminating the sign problem, and the calculation of 
2fc(T,N) is not easy actually. We rather use a density of state method with fixing the plaquette 
variable and apply an approximation proposed in [Q] to avoid the sign problem. 

We introduce a probability distribution function of the plaquette, which is defined by 

W(P',P,li q ) = J mi 8(P'-P) (detM)% 6 ^ p , (3.1) 

where 8 (x) is the delta function. For later discussions, we define the average plaquette P as P = 
— 5 , g/(6j8A^si te ) and the quark matrix M as independent of j8. A^ s i te is the number of sites. The 
plaquette distribution functions for \i q = are shown in Fig. [j] (right) for each j8. We denote 
w(P,p) = W(P, j8,0). The normalized partition function is rewritten as 

^ffi = Jw(P,P,Li q )dP = jR(P,ti q )w(P,f3)dP. (3.2) 

Here, R(P,jj, q ) is the reweighting factor for finite fj, q defined by 

R(p , „ v _ J2U S(P'-P)(detM( N )) N t n ggojj; 1^=0) 

1 f@US(P>-P)(d e tM(0)) N < (S(P'-P)) { ^= ) ' 1 j 
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This R(P,n q ) is independent of /3, and R(P^\x q ) can be measured at any /3. In this method, all 
simulations are performed at \i q = and the effect of finite \i q is introduced though the operator 
detM(jU 9 )/detM(0) measured on the configurations generated by the simulations at \i q = 0. 

Because QCD has time-reflection symmetry, the partition function is invariant under a change 
from jj, q to — }i q , i.e. R(P, —jl q ) = R(P, fl q ). Moreover, the quark determinant satisfies detM(— fi q ) = 
(detM(jU*))*. From these equations, we get 

[R(P,ll q )]*=R(P,Ll* q ). (3.4) 

This indicates that R(P,n q ) is real in the case of real \i q , i.e. \i q = \i* q . Then, the probability 
distribution of the plaquette given by R(P,n q )w(P,fi) is real. However, once the imaginary part of 
H q becomes nonzero, R(P,n q ) is not a real number any more. We thus write the partition function, 



Je i ^\R(P^ q )\w(P,l5)dP= fe^ p ^\R(P^ q )\w(P,f5)(^j rff (3.5) 



Because this complex phase (p vanishes at = 0, (f) does not have the -independent part, and 
the problem of the statistical error due to the /ii m -independent part is eliminated. As we discussed 
in the previous section, if \R(P,jj, q )\w(P,p)(d(j) /dP)~ l becomes a double-peaked function of (j) and 
the distance between these peaks is equal to (2n+ 1)71 with an integer n, the partition function 
becomes zero, i.e. a Lee- Yang zero appears. To investigate whether the probability distribution has 
double-peak or not, we introduce an effective potential defined by 

V eff (f = -hi\R(P^),n q )\-]nw(P((l>),P)+]n(^) (P(<j>),H q ). (3.6) 

\ a " / norm 

Here, we normalize the value of d<p/dP by the maximum value for each \i q . This effective potential 
is a function of P in practice. We define Ves(P, j3 , pL q ) = V e ff (0 (P) , j3 , /^^ ) . 

Performing Monte-Carlo simulations, we calculate these three quantities, \R(P,jx q )\, (j>(P,[i q ) 
and w(P,P). However, the exact calculation of the quark determinant is difficult except on small 
lattices. In this study, we estimate the quark determinant from the data of Taylor expansion coeffi- 
cients up to 0(n q ) around jj, q = obtained by a simulation of two-flavor QCD with p4-improved 
staggered quarks in [J3j] . The truncation error has been discussed in [Q]. Because we define 6 by 
the Taylor expansion of IndetM, 6 is not restricted to the range from —% to 7t. 

Next, we discuss the sign problem in the calculation of V e ff. We denote the quark determinant 
as iVfln[detM(jU 9 )/detM(0)] = F + id. Histograms of 8 are shown in Fig. |l] (left). We fitted the 
histograms to Gaussian functions. The results are the dashed curves. The distributions seem to be 
well-approximated by Gaussian functions. Here, we perform a cumulant expansion, 



(exp(F + i0)) = (e F )exp 



i(fl) - \mf) -^((Ad?) + i(AFA0) - ^(AF (Ad) 2 ) + • 



(3.7) 



with AX = X — (X). If the distribution of 6 is of Gaussian, the 0(6") terms vanish for n > 2 in 
this equation [Q]. Moreover, since d ~ 0(n q ) and F ~ 0(\l\), this expansion can be regarded as 
a power expansion in \i q . If the expansion of Eq. ( [377| ) is good, we can extract the phase factor 
e lif from R easily and the sign problem in \R\ is eliminated. We deal with the first two terms, i.e. 
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Figure 2: The derivative of the effective potential dVes/dP at j3 = 3.65 for (pL lm /T) 2 = 0.0 (left), 
(/W r ) 2 ~ (fiim/T) 2 = 1.0 (middle) and 4.0 (right). We measured them at the peaks of the plaquette 
histograms in Fig. [j] (right) and interpolated the data by a cubic spline method. 



i(6) and — ((Ad) 2 } /2, assuming the Gaussian distribution. The correlation terms between F and 6 
are also neglected as a first step. Because we calculate the expectation value with fixed P and the 
values of F and P are strongly correlated, the AF may be small once P is fixed. Then, » (8). 

The results of the first derivative of V e ff(P, jS,/!^) are shown in Fig. f| for various H q /T with 
(Mim/r) 2 = (left) and with Re[(^l q /T) 2 } = (^i Re /T) 2 - (m m /T) 2 = 1 (middle) and 4 (right). 
j8 = 3.65 is adopted for these results. We discuss dV e s/dP instead of V e & itself because d 2 V e s/dP 2 
is independent of j3. dV^jdP at different j8 can be estimated by the equation, 

^rttfl) = ^rttAO -6(J8 - jSoKite, (3.8) 

under the parameter change from /3o to /3. This equation is derived from the equation, w(P,fi) = 
e 6 W-W NsiaP w(P,fio), which is given from the definition, Eq. (j3?T|). Therefore, the form of dV eS /dP 
as a function of P does not change with /3 up to a /3 -dependent constant. Using the behaviors of 
dV e s/dP, the value of P minimizing V e ff, i.e. dV e ff/dP = 0, can be also controlled by /3. 

If the effective potential Veff is a double- well function of <j>, dV e s/d<f> is an S -shaped func- 
tion and vanishes three times. In such a case, there exists a region of P where the derivative 
of dV e f[/dP = (dV e ff/d<p)(d<p/dP) is negative, since <p is a monotonically increasing function of 
P. The left panel of Fig. ^| shows that the region of d 2 V e fi/dP 2 < appears at high density, i.e. 
(/iRe/r) 2 > 6 for jUim/T = 0. Also, in the region of large jX\ m /T, dV^/dP becomes an S-shaped 
function, which is shown in the middle and right panels of Fig. § for (/iR e /r) 2 = 1 and 4. 

Next, we investigate the boundary at which dV e ft/dP changes to an S-shaped function from a 
monotonic function, which is shown as a line in Fig. |3| (left). On this line, d 2 V e ff/dP 2 = at a value 
of P. The constant part of dV^/dP is changed by /3. The values of /3 at which both d 2 V e n/dP 2 and 
dV e g/dP vanish simultaneously are indicated in this figure. Above this line (large Im(/i 2 ) or large 
Re(/i 2 )) with this /3, the double-well potential appears, and the phase cancelation occurs near this 
line. The temperatures corresponding to these values of j8 are shown above the curves of w(P,f3) 
in Fig. [T] (right). This result suggests the existence of singularities (Lee- Yang zeros) in the region 
of large Im(/x 2 ) as well as the region of large Re(/i 2 ), and the boundary at large Im(/x 2 ) is closer 
to ji q = than that at large Re(/i 2 ). The distance to the boundary from \i q = is shown in Fig. || 
(right) for each j8 . The blue circles are the distance to the line in the complex plane and the red 
square is that on the real axis. The boundary on the real axis is essentially the same as the critical 
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Figure 3: Left: Boundary where V e g changes to double-well type in the complex (jiq/T) 2 plane. Below 
this line, V e s is always of single-well. The numbers in this figure are the values of j3 at which dV e s/dP 
and d 2 V B ff/dP 2 vanish simultaneously. Right: The distance to the boundary form ji q = for each /3. This 
corresponds to the radius of convergence. The red symbol is the results on the real axis. 

point in This distance to the boundary is approximately equal to the convergence radius of an 
expansion of ln^/(Vr 3 ) = Y M C2n{V- q /T) 2n . 

4. Conclusions 

We studied the singularities of QCD in the complex \i q plane by a numerical simulation. Be- 
cause the sign problem makes the calculation of the partition function difficult, we discussed a 
probability distribution function of a complex phase instead of the partition function itself. In this 
calculation, we used a kind of reweighting method together with the approximation in [Q]: the 
quark determinant is estimated by a Taylor expansion, and the Gaussian distribution of the com- 
plex phase of detM is assumed. We found that there is a region where the phase distribution has 
two peaks, suggesting the existence of singularities, in the region of large Im(/i^) as well as of large 
Re(jU^). We moreover estimated the distance to the nearest singularity from \i q = in terms of pL q 
for each temperature. The distance is regarded as the convergence radius of a Taylor expansion of 
the thermodynamic potential. Although our simulation is performed with relatively heavy quark 
mass, the result suggests that convergence radius at the temperature of the real QCD critical point 
(T cp ) may be longer than the convergence radius in the crossover region at T > T cp . 
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